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We calculate the renormalized effective two-, three-, and four-body interactions for N neutral 
ultracold bosons in the ground state of an isotropic harmonic trap, assuming two-body interactions 
modeled with the combination of a zero-range and energy-dependent pseudopotential. We work 
to third-order in the scattering length a t (0) defined at zero collision energy, which is necessary to 
obtain both the leading-order effective four-body interaction and consistently include finite-range 
corrections for realistic two-body interactions. The leading-order, effective three- and four-body 
interaction energies are U 3 (u) = -(0.85576. ..)[a t (0)/a(uj)} 2 + 2.7921(1) [a t (0)/a(uj)} 3 + O(at) and 
U 4 (lu) = +(2.43317.. .)[a t (0)/cr(V)] 3 + 0(a$), where u and a(u) are the harmonic oscillator fre- 
quency and length, respectively, and energies are in units of fvuj. The one-standard deviation error 
±0.0001 for the third-order coefficient in Uz(ui) is due to numerical uncertainty in estimating a 
slowly converging sum; the other two coefficients are either analytically or numerically exact. The 
effective three- and four-body interactions can play an important role in the dynamics of tightly con- 
fined and strongly correlated systems. We also performed numerical simulations for a finite-range 
boson-boson potential, and it was comparison to the zero-range predictions which revealed that 
finite-range effects must be taken into account for a realistic third-order treatment. In particular, 
we show that the energy-dependent pseudopotential accurately captures, through third order, the 
finite-range physics, and in combination with the multi-body effective interactions gives excellent 
agreement with the numerical simulations, validating our theoretical analysis and predictions. 

PACS numbers: 31.15.ac,31.15.xp,05.30.Jp,67.85.-d 



I. INTRODUCTION 

Effective multi-body interactions arise when quantum 
fluctuations dress the intrinsic interactions between par- 
ticles. They play a central role in quantum field theories 
and exemplify the significant difference between inter- 
actions in classical and quantum theories. For example, 
even for a quantum field that has only intrinsic two-body 
interactions at high energies, at low-energy scales, after 
the high-energy degrees of freedom are coarse-grained 
away, the field will manifest at some level effective n- 
body interactions. The ability to trap and control sys- 
tems of ultracold neutral atoms [TJ [2] has created new 
opportunities to study this physics in the laboratory. Ef- 
fective three-body interactions in the limit of large two- 
body scattering length have in particular received a great 
deal of attention, motivated both by the predictions of 
universal behaviors [3H2] and the ability to use ultracold 
atoms to study physics ranging from molecular [10] to nu- 
clear scales [TTJ [12] . Recently, attention has focused on 
Efimov-like states and universal behaviors for four-body 
systems, again in the limit of large scattering lengths 

[IMS]. 

Here, we focus on the opposite regime of weakly in- 
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teracting neutral bosons with small scattering lengths. 
Even in this limit, effective higher-body interactions can 
be important, particularly for tightly confined or strongly 
correlated particles. This is seen dramatically in [T7] , 
where a superfluid of bosonic atoms is quenched by sud- 
denly increasing the depth of an optical lattice. After the 
quench, which creates a non-equilibrium state of strongly 
correlated bosons, beating effects due to multiple distinct 
interaction energies, as expected from effective three- and 
higher-body interactions p~8j [19] , are seen in the collapse 
and revival oscillations of the first-order coherence. Ef- 
fective multi-body interactions should also have played 
a role in previous collapse and revival experiments [20l - 
[22] . although in those cases inhomogeneities may have 
masked their signature. More recently, effective three- 
and four-body interactions have been used to demon- 
strate atom-number sensitive photon-assisted tunneling 
in optical lattices [23], and their influence has been seen 
in precision measurements on Mott-insulator states of ul- 
tracold atoms [24] . A number of studies also suggest that 
elastic multi-body interactions can play an interesting 
role in generating exotic quantum phases in optical lat- 
tices or modifying the superfluid to Mott-insulator phase 
transition [25H32] . 

In this paper, we use renormalized quantum field the- 
ory [33] to calculate the perturbative ground-state energy 
for N ultracold neutral bosons in a three dimensional 
isotropic harmonic potential with angular frequency cj, 
and extract from it the effective m-body interaction en- 



ergies ^(cj), Us(lj), and U±{uS) as a function of u. The 
key purpose of the present paper is to (i) systematically 
develop a renormalized quantum field theory approach 
for ultracold trapped bosons including finite-range ef- 
fects, (ii) determine the leading-order four-body interac- 
tion, and (iii) validate the formalism through comparison 
with numerical results. To obtain effective four-body in- 
teraction energies it is necessary to work through third 
order in the two-body scattering length. We use renor- 
malized perturbation theory (see [33J), which develops 
an expansion around physical as opposed to bare cou- 
pling parameters, to systematically cancel the multiple 
divergences that arise at higher-orders in quantum field 
perturbation theory. (In this paper, the physical coupling 
parameter is defined in terms of the measured scattering 
length, or alternatively the measured energy shift, for two 
interacting ultracold boson in a harmonic trap at a spec- 
ified trap frequency.) Renormalized perturbation theory, 
which is more commonly used in high-energy physics, in 
this context naturally describes how the effective inter- 
actions depend on trap frequency. An example of the 
power of renormalized perturbation theory to capture 
low-energy physics is that we independently reproduce, 
through third order, the two-body ground-state energies 
calculated in [34 . More fundamentally, the analysis in 
this paper provides an explicit example of renormaliza- 
tion physics and running coupling constants that can be 
directly probed using trapped ultracold bosonic atoms, 
and used to test central concepts in effective field theory. 

To calculate effective interactions for confined bosons, 
we first assumed that the two-body interactions could be 
described in the low-energy, s-wave limit by an energy- 
independent zero-range ^-function pseudopotential. To 
test our perturbative predictions, we then numerically 
calculated TV-boson ground-state energies using a finite- 
range two-body Gaussian model potential. Comparison 
with the numerical results revealed that finite-range ef- 
fects must also be taken into account for an accurate 
description of realistically interacting bosons. In this pa- 
per, we show that both the finite range effects and ef- 
fective interactions are accurately captured by the com- 
bination of zero-range and energy- dependent ^-function 
pseudopotentials. Including the finite-range corrections, 
we are able to validate our analytic and numerical cal- 
culations of all perturbation theory coefficients through 
third order. 

The basic idea in our approach is the following: we "in- 
tegrate out" excited vibrational states thereby trading a 
multi-orbital theory with intrinsic two-body interactions 
for a single-orbital theory with effective multi-body inter- 
actions. The latter can provide a simple but powerful al- 
ternative description of the low-energy few-body physics. 
The quantum fluctuations to excited states both dress the 
two-body interactions and generate effective higher-body 
interactions. The idea is illustrated in Fig. [I] We showed 
in [HI [19] how this approach can be used to approxi- 
mately incorporate the influence of higher bands via the 
simple modification of adding higher-body interactions 
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FIG. 1: Illustration of the idea of replacing a multiple-orbital 
(or multiple-band) model with orbitals i = 0, 1, ... with only 
intrinsic two-body interactions by a single-orbital (or single- 
band) model with effective multi-body interactions between 
renormalized states. A ground state multi-body model can 
be useful when virtual excitations of bosons to excited vibra- 
tional levels are important. 



to the single-band Bose-Hubbard model [35j [36] . 

Beyond applying directly to ultracold neutral bosons 
in an isotropic harmonic potential, our results can give 
qualitative insight into the effective interactions for other 
trapping potentials. They can also be used for rough 
approximations to the effective two-, three-, and four- 
body interactions in anisotropic potentials, and for neu- 
tral bosons in optical lattices. In the latter case, however, 
anharmonicities are important. For example, we esti- 
mate an approximately 30% anharmonic correction to the 
three-body interactions for 87 Rb in typical lattices. The 
role of anharmonicities for collapse-and-revival dynamics 
in optical lattice systems has been analyzed further in 
[37] . Inhomogeneities and the effect of a background har- 
monic potential on lattice collapse-and-revival dynamics 
has been studied in J38[ [39] . 

Tunneling also has an influence on collapse and re- 
vival in optical lattices [40-42 . In deep (post-quench) 
lattices the typical tunneling energy is nearly an order 
of magnitude smaller than the effective three-body inter- 
action energy, making the latter effect dominant. Tun- 
neling should, however, be of comparable importance to 
the effective four-body interactions. Approaches apply- 
ing effective interaction methods to tunneling in lattice 
or multi-well systems include [43-46], and related meth- 
ods for analyzing physics involving interactions, corre- 
lations, higher bands, and quantum tunneling in lattice 
systems include [47- 52 . Fermionic systems and fermion- 
boson mixtures also yield interesting types of effective 
interactions that have received increasing attention (e.g., 
[3T| [37j I53H56] ). as well as three-body interactions of 
fermions and polar molecules in lattices [25] . 

For experiments with 87 Rb at typical lattice densi- 
ties the recombination rate [5j [57] is one or more or- 
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ders of magnitude smaller than the frequencies associ- 
ated with both the effective three- and four-body ener- 
gies, and therefore the elastic effective interactions de- 
scribed in the present paper are more important than in- 
elastic multi-body interactions driving loss. Roughly, we 
expect three-body recombination to scale at fourth order 
in the scattering length [58] , and in the future we would 
like to understand both elastic and inelastic interactions 
in a unified framework. The role of effective three-body 
interactions in thermalizing a homogenous ID Bose gas 
has also been studied [59], and it would be interesting 
to investigate this physics in the context of a 3D optical 
lattice system. 

The remainder of this paper is organized as follows. In 
Sec. [IlJ we provide an overview of our results. Section III 
compares the perturbation theory predictions to numer- 
ical estimates for finite-range interactions. Sections IV 
and |V| describe the details of the renormalized perturba- 
tion theory used to obtain the effective multi-body inter- 
actions. Section HVl defines the renormalized Hamiltonian 
and derives the first- and second-order corrections, while 
Sec. [V] derives the two-, three-, and four-body interac- 
tion energies through third order. Section VI summarizes 



our results and conclusions. Finally, the appendices give 
derivations of a number of technical results used in the 
paper. 



II. OVERVIEW 

We find the effective interactions of N ultracold bosons 
in the ground state of an isotropic harmonic oscillator 
with pairwise interactions modeled by a zero-range 5- 
function pseudopotential 



.) = 52< 5(3) (r ._ r . ); 



(1) 



where is the position vector of the i th boson. We 
assume there are no intrinsic three- or higher-body in- 
teractions. The two-body coupling constant #2 is re- 
lated to a t (0), at first order in perturbation theory, by 
g 2 = 4tt (h 2 /rriA) a t (0) + (9([a t (0)] 2 ), where tua is the bo- 
son mass, a t (0) is the physical s-wave scattering length 
measured in the limit that the trap frequency and colli- 
sion energy go to zero, and O([a t (0)] 2 ) are terms of order 
[a t (0)] 2 and higher. At higher orders, the relationship 
between #2 and a t (0) is modified, and in Sees. IV and 
|V| we generalize the perturbation theory as an expan- 
sion around the physical trap scattering length a t (uo) 
defined for a harmonic potential with frequency ujq. In 
this overview, we summarize our results to third order in 
a t (0), i.e., the special case uq = 0. 

We obtain the ground-state energy of N bosons as 
an expansion E = En=o^ (n) > where eM 

is propor- 
tional to [at(0)] n - Throughout this paper energies are ex- 
pressed in units of the harmonic oscillator energy huo. The 
zeroth-order (one-body) energy is E^\uj) = s N, where 
eo = 3/2 is the dimensionless single-particle ground-state 



energy. The n th -order energies for n > can be expanded 
as 



N 



(2) 



where ( m ) is the binomial coefficient. The sum goes up 
to the minimum of N and n + 1, and the n th -order con- 
tributions to the m-body interaction energies (in units of 
fyjo) are 



U M (cj)=c M ( at(0)V 



(3) 



where the harmonic oscillator length for an isotropic po- 
tential with frequency u is 



Table [i] gives the values of cfiP obtained in Sees. 



(4) 



(3) 



IV 



and |v[ The two-body coefficients cp , cp , and 
dependently reproduce the results in [34 , if the exact 
solution found there is expanded through third order. 
The coefficient c\ , in particular, is nontrivial and pro- 
vides a strong consistency check that the renormalized 
perturbation theory captures the two-body low-energy 
interactions correctly. 

(2) 

The analytic value of the three-body coefficient c\ was 

previously found in [18] . The coefficient found here 
extends that result to third order in a t (0). The value of 

cP given in Table [i] combines both analytic and approxi- 
mate numerical results, and the uncertainty is due to the 
slow convergence of one of the numerically determined 
sums (see App. B2). 

(s) 

We also obtain the coefficient c\ , which gives the 

leading-order contribution to the effective four-body en- 
ergy. The coefficient cP combines numerical and ana- 

lytic results, but unlike c% ) has high precision because of 
the fast convergence of all the contributing terms. Note 
that cP and cP have similar magnitudes, and conse- 
quently we need to include the effective three-body cor- 
rections when effective four-body effects are important 
or of interest. At the end of Sec. [II] we show that the 
correction from the third-order terms becomes significant 
for ultracold atoms in trap potentials with relatively tight 



confinement. The coefficients cP and c^ J have not pre- 
viously been reported in the literature. 



,(3) 



In Sec. Ill we compare the predictions for zero-range 
interactions to numerical calculations for a Gaussian 
boson-boson interaction potential and find significant ef- 
fects from its finite-range nature. We show that these are 
accurately modeled by adding to the zero-range pseu- 
dopotential V2 an energy- dependent (higher-derivative) 
pseudopotential [11] 



-f[^/ (3) (r 



1 )+S^(r i -r J )^l], (5) 
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(2/?r)(l - log 2) = +0.19535.. 



Effective Interaction Energy Coefficients 
Two- 

4 1} = (2/tt) 1 / 2 = +0.79788.. 

c (2) 

c^ 3) = (2/tt) 3 / 2 (1 - -3 log 2+ flog 2 2) = -0.39112... 
4 1,2) = (3/4)(2/tt) 1 /2 = +0.59841... 

4 2) = (2/7r){-4v / 3 + 6[l-21og2-log(2- v^)]} 
= -0.85576... 

4 3) = -12(2/tt) 1 / 2 (1 - log2)4 2) + 12c4 3) - 6a^ - 18a^ 3) 
= +2.7921(1) 

Four-body 

4 3) = 48c4 3 l + 48a^ - 72<4 3) = +2.43317... 

TABLE I: The coefficients and d^ ,2 \ which give 

the n th -order correction to the ra-body effective interac- 
tion energies Um^ (w) [see Eqs. ^ and ^] for neutral 
bosons in an isotropic harmonic potential. The results for 
c£\c£\c£\df£ ,2 \ and are exact. The coefficients 4^ 
and cf^ are given in terms of parameters a^\a^\ etc., de- 
fined in Table [n) We have obtained exact analytic expressions 
for ot^ , af\ , and af^ . The numerical approximations for af\ 

(3) 

and ay are obtained to very high precision, but slow con- 
vergence of the expression giving is responsible for the 



(3) 



uncertainty in the value of c. 



which has been symmetrized to make it Hermitian. The 
operators V;j and are gradients with respect to the 
relative separation — r^, acting to the left and right, 
respectively. The coupling constant is 

g' 2 = (to—) ( lr eS [at(0)} 2 ) + OOeffMO)] 3 ), (6) 



mA / \2 

where r e ff is the effective range jgO] • To first-order in g' 2 , 
the shift to the iV-body ground-state energy is 



£(1,2) 



with 



(1,2) / r e ff 



flt(O) 



(7) 



(8) 



The superscript (1,2) indicates that the term is first order 
in r e ff and second order in a t (0), and 4 is given in 
Table |H 

The potential V2 is proportional to r e ff [a t (0)] 2 /<j(o;) 3 

and we consider in this paper a regime where a t (0) ~ 

fi 2) 

r e ff <C ct(cj), such that V 2 and therefore C/g ' (^) can be 
treated as if the contribution is third order in a t (0) . This 
approach is supported by the comparison between the 
perturbative energies and the energies for the Gaussian 



potential with spatial widths Tq < 0.01cr(o;) in Sec. Ill 



fi 2) \ 

Adding the contribution £7 2 ' i 00 ) to ^ ne two-body inter- 
action energy extends our results to more realistic sys- 
tems, like ultracold atoms that interact through finite- 
range van der Waals potentials. 



Equation ^ organizes the TV-body energy in powers 
of the free-space s- wave scattering length a t (0). Alter- 
natively, combining our results, we can reorganize the 
energy in terms of m-body contributions as 



E = eoN+-U 2 (u)N(N-l) 

+ ^U 3 (uj)N(N-l)(N-2) 



(9) 



+ -U 4 (uj)N(N - 1)(N - 2)(N - 3) + 

where through third order the two-body interaction 
energy is 



(!) /Ot(0)\ (2) f°t(0)\\(3) ^t(0) x3 



(j{uj) 



,(i,2) ( r^_\ ( q t (0) 



O 



[at(0)} 4 



o 



the three-body interaction energy is 



(10) 



a(ui) J ' ~ 3 \a(ui) 
r e ff[a t (0)] 3 



/Mf 



O , 

and the four-body interaction energy is 



(ii) 



U 4 (lo) = c 
+ 



(3) /Ot(0) 



4 

r e ff[a t (0)] 3 

[<T( W )]4 



O 



[«t(0)] z 



(12) 



The four-body interaction energy [^(cj), although 
comparatively small, can lead to qualitatively impor- 
tant effects, particularly for traps with stronger con- 
finement. For example, for N = 4 87 Rb atoms and 
a t (0)/cr(uj) = 0.05, corresponding to a 10 4 Hz trap fre- 
quency, the four-body energy should generate a distinct 
approximately 60 Hz beating frequency in collapse-and- 
revival oscillations, using our harmonic trap results to 
estimate the energy in an optical lattice potential. These 
effects should be measurable as long as tunneling and 
trap inhomogeneities are sufficiently reduced [T7] . 



Using the effective interaction energies in Eqs. (10), 



(11), and (12), we can construct a single-orbital effective 



Hamiltonian 



^— ' ml 

m=2 



(13) 



where a (aJ) annihilates (creates) a boson in a renormal- 
ized single-particle ground state. The effective Hamil- 
tonian can be used to incorporate some higher-band 
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physics, via effective multi-body interactions, into a 
single-band Bose-Hubbard model [T8] . 

The effective interaction energies can be tuned by 
changing either the scattering length a t (0), for example 
with a Feshbach resonance [2 , or the trap frequency uj of 
the confinement [6lJ[62]. For example, for a fixed a t (0), 
this tuning follows from rewriting the U m (uj) in terms of 
the characteristic scattering energy ftw s = h 2 / mp\at(^S)} 2 . 
That is, we write U m (uj) = U m (uj)(uj/uj s ) such that 



20 kHz 40 kHz 60 kHz 



i 1 ) 



(uj/uj s ) 



3/2 , J2) 



(uj/uj s ) 2 



(14) 



+ (4 3) + 4 X ' 2) [r eff /a t (0)])(^M) 5 / 2 + 0[(uj/uj s )% 
U 3 (uj) = cf\u/u s ) 2 + c^(uj/uj s )^ 2 + 0[(uj/uj s )% 

(15) 

U 4 (uj) = cf(uj/uj^ 2 + 0[(uj/uj s ) 3 }. (16) 

Figure [2] shows, for the case of a zero-range poten- 
tial (i.e. r e ff = 0), the two-body energies U^fa) (in 

~ (2) ~ (3) 

the inset) and (cj) + ^ (cj), the three-body ener- 

~ (2) ~ (2) ~ (3) 

gies C/3 (uj) and C/3 ; (cj) + U% (uj), and the four-body 



energy [/"^(u;) versus uj/uj s . As expected, L 7 ^ 1 " 1 ^) is the 

~ (2) 

largest contribution. The line labeled U% (uj) shows the 

second-order three-body result found previously in [18], 

(2) ~ (2) ~ (3) 

due to the c 3 ; coefficient, and the line (uj) + /7 2 (uj) 

shows the scale of the correction from the third-order co- 
efficient C3 . The effective three- and four-body energies 
have opposite signs and are of similar magnitude. Fi- 
nally, the line labeled /7f xact (cc;) — U^\uj) shows the good 
agreement with the exact two-body results from [34] for 
the regularized zero-range potential. 

It is interesting to directly compare the relative sizes 
of the second- and third-order corrections for 87 Rb 
in a trap. For small magnetic field strengths, the 
87 Rb scattering length and effective range are approx- 
imately 5.3 nm and 7.9 nm, respectively [2J. For a 
trap frequency of 10 2 Hz, and thus a t (0)/a(uj) = 0.005 
("weak" confinement), the third-order two-body terms 
4 3 Vt(0)/a( W )] 3 and4 1 ' 2) [a t (0)/cr( W )] 2 [r eff (0)/a( W )] are 
1% and 2% of the second-order two-body contribution 
c>P [a t (0) /a(uj)] 2 . Similarly, the third-order three- and 

four-body terms 4 3) M0)/cr(uj)} 3 and cf } [a t (0)/a(u)} 3 
are each about 1.5% of the second-order three-body con- 
tribution cf\a t (0)/(j(uj)] 2 . 

For a trap frequency of 10 4 Hz, and thus a t (0)/a(u) = 
0.05 ("strong" confinement), the third-order two-body 
terms increase giving approximately 10% and 20% cor- 
rections compared to the second-order two-body con- 
tribution. Similarly, the third-order three- and four- 
body terms increase giving approximately 15% correc- 
tions compared to the second-order three-body contri- 
bution. (Notice, however, that the third-order effective 
two-body coefficient and the finite-range coefficient have 
opposite signs, and hence their contributions partially 
cancel.) In typical optical lattice collapse- and-revival ex- 
periments with 87 Rb the confinement is even stronger 
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FIG. 2: (Color Online) Perturbative predictions for dimen- 
sionless effective interaction energies U m (oj) versus oj/uj s for 
fixed scattering length at(0), in units of the energy huj s = 
h 2 /mA[at(0)] 2 - The inset shows the first-order two-body en- 
ergy tj^(uo). The main figure shows the second- and third- 
order corrections to the two-, three-, and four-body energies, 
assuming no finite-range corrections. The top and right axes 
in both figures show the energies converted to frequency units 
by multiplying by u s /2tt, assuming 87 Rb with a t (0) = 5.3 
nm, m A = 86.9 u, and uj s /2ty = 4.14 MHz. The line labeled 
L/| xact (w) — TJ^ (uj) gives values using the exact two-body re- 
sults for U2(uj) from |34| . 



and the ratio a t (0)/a(uj) is on the order of 0.05 — 0.10 
[17, 20-22 . In this regime we expect non-perturbative 
effects to also become increasingly important. 



III. COMPARISON OF PERTURBATIVE 
ENERGIES WITH ENERGIES FOR 
FINITE-RANGE INTERACTIONS 

This section compares the predictions of the pertur- 
bative ground-state energies for a zero-range ^-function 
interaction potential, summarized in Sec. [XT] and de- 
rived in Sees. [IV] and |V| and numerically obtained en- 
ergies for TV-boson systems with finite-range interac- 
tions. We show that the leading-order contribution 
of an energy-dependent pseudopotential accurately cap- 
tures the finite-range effects, and allows us to also vali- 
date the analytic and numerical coefficients found from 
the zero-range perturbation theory. 

We use a finite-range interaction model based on a 
Gaussian two-body potential V g (r) = Vo exp[— (r/ro) 2 /2] 
with depth (or height) Vq and width ro [63j[64]. For a 
given width ro, we adjust the depth Vo such that V g (r) 
produces the physical free-space s- wave scattering length 
a t (0) at zero collision energy. We restrict ourselves to 
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FIG. 3: (Color online) Free-space scattering quantities for 
the Gaussian model potential, with all lengths expressed 
in units of ct(cj). Circles and squares show the volume 
^eff[at(0)] 2 /[cr(a;)] 3 (in the main figure) and the effective range 
r e s/<j(oj) (in the inset) as a function of a t (0)/a(uj) for the 
Gaussian potential with ro = 0.005ct(cj) and ro = O.OIct(cj), 
respectively. 



depths Vo for which V g supports no two-body s-wave 
bound state in free-space. This implies that Vo is positive 
for a t (0) > and negative for a t (0) < 0. 

An energy-dependent free-space scattering length for 
two particles with relative energy E re \ and relative wave 
number k re \ = ^/mxE re \/h can be defined as 

tan((Sf(fe re i)) 



^f(^rel) = 



k 



rel 



(17) 



where 5f(fc re i) is the free-space s- wave phase shift. The 
effect of a finite-range potential on the free-space scatter- 
ing of two ultracold bosons can be captured by Taylor- 
expanding 5f(fc r ei) [60, 65 , giving 



at(E rel ) = Ot(0) + -r eff [a t (0)] 2 /c : 



rel 



(18) 



where r e R is the effective range parameter which describes 
the lowest-order energy dependence of the phase shift 

HUES]. 

Figure [Hi] shows the effective range r e ff and the 
"volume" r e ff[a t (0)] 2 for two bosons interacting with 
the Gaussian potential with two different choices of 
ro/a(u) <C 1. (The volume factor here characterizes the 
leading-order effective-range correction to s- wave scatter- 
ing.) We extract r e ff by fitting the numerically evaluated 
— tan(5f(fc re i))/fc re i to the right-hand-side of Eq. ( p~8| ) for 
small scattering energies. The effective range is positive 
for negative a t (0), negative for small positive a t (0), and 
diverges as a t (0) — » 0. Importantly, since a t (0) = 
implies Vo = (no scattering potential), the volume 
r e ff[a t (0)] 2 also vanishes when a t (0) = 0, as seen in the 
main part of Fig. Ill The divergent behavior of the ef- 
fective range is also observed for realistic van der Waals 
potentials [66] and indeed for any potential that falls off 



faster than 1/r 5 [65 , although for these potentials (un- 
like the Gaussian) r e ff[a t (0)] 2 is finite but non-zero in the 
limit a t (0) -> 0. 

We determine the ground-state energy of N = 3 and 
N = 4 bosons interacting through the Gaussian model 
potential under external spherically symmetric harmonic 
confinement using a basis set expansion that expresses 
the relative iV-body wave function in terms of explicitly 
correlated Gaussians 64 



N h 



Iprel = ^2 UkS eXP 



k=l 




(19) 



The Uk denote expansion coefficients, A5 is the number 
of basis functions, and S symmetrizes the wave func- 
tion under the exchange of any pair of bosons. The 
Ni>xN(N—i)/2 variational widths v\^\ chosen stochasti- 
cally from the interval [ro/5, 4cr(u;)], are optimized semi- 
stochastically following the scheme outlined in Ref. [64] . 
In brief, the variational method works as follows. Assume 
we have a basis set consisting of j — 1 basis functions that 
yields a ground-state energy estimate Ej-±. To add the 
j th basis function (j < A/5), we generate a few thou- 
sand trial functions. For each trial function, we solve 
for a trial ground-state energy by diagonalizing a j x j 
dimensional generalized eigenvalue problem. (It is a gen- 
eralized eigenvalue problem because the basis functions 
are nonorthogonal.) We choose as the j th basis function 
the one which makes Ej smallest, and repeat this process 
for the (j + l) th basis function until j = A5. A key ben- 
efit of the explicitly correlated basis functions is that the 
Hamiltonian and overlap matrix elements have compact 
analytical expressions [64] . 

Convergence is analyzed by investigating the depen- 
dence of the energies on Nb and by performing calcula- 

(k) 

tions for different sets of widths v\- . To meaningfully 
compare numerical three- and four-body energies i?FR 
for the finite-range (FR) interaction potential with per- 
turbative results up to order [a t (0)] 3 , the numerical accu- 
racy of the finite-range energies should be notably better 
than |a t (0)/cr(cj)| 3 . For example, for |a t (0)| = 0.001a(cj) 
and I at (0)| = 0.01a(a;), this implies numerical accuracy 
better than 10 -9 and 10 -6 , respectively. An analysis of 
the basis set error shows that our Af-body energies are 
sufficiently accurate to test the perturbative predictions 
up to order [«t(0)] 3 for |a t (0)| > O.lro, using about 100 
and 500 basis functions for N = 3 and N = 4, respec- 
tively. Our numerical accuracy is insufficient to test the 
perturbative predictions for smaller |a t (0)|. 

Figure E] shows the quantity [E FR - E^ - E&] x 10 4 
versus a t (0)/cr(cj), with the finite-range energies E^r 
numerically computed using ro = 0.01<t(co>) (the blue 
squares) and ro = 0.005<r(co>) (the red circles). We have 
subtracted the energies E^ and E^ obtained from the 
perturbative theory to better examine the physics beyond 
first order in a t (0). The solid line is [E^ + £ (3) ] x 10 4 
from the perturbative theory with r e ff = 0. Panels (a) 
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0.01 



at(0)/a(co) 



FIG. 4: (Color online) The quantity [E FR - E w - E (1) ] x 10 4 
versus at(0)/cr(cj) for N = 3 and 4 in panel (a) and (b), 
respectively. (Energies are in units of hw.) The finite-range 
energies Efr are numerically computed with ro = 0.01cr(ct;) 
(blue squares) and ro = 0.005ct(cj) (red circles). The solid 
line is [E^ + E^] x 10 4 found from the perturbative theory 
with the zero-range potential. 



and (b) give the energies for N = 3 and 4 bosons, respec- 
tively. For N = 3, we see that finite-range corrections 
to the zero-range theory become more significant for in- 
creasing ro- 

In Figs.[5ja) and (b), we multiply the N = 3 and 4 en- 
ergies E FR - - by [a(u)/a t (0)} 2 . The perturba- 
tive predictions for (E^ +E^)[a(uj) / a t (ti)] 2 are straight 
lines. The nonperturbative numerical results are for po- 
tentials with ro = 0.005<t(o;) and ro = O.Ola(cj). The fig- 
ures show that the scaled numerical results are singular 
near zero scattering length, and only approach the zero- 
range perturbative results with increasing |a t (0)|. More- 
over, by decreasing ro the difference between the pertur- 
bative results and the scaled finite-range energies is re- 
duced, and we conclude that the divergences at a t (0) = 
are due to the finite range of the Gaussian potential. Mul- 
tiplying the energies by [a(u)/a t (0)] 2 has magnified the 
finite-range corrections, showing that an effective field 
theory description for finite-range potentials requires cor- 
rections to the zero-range (5-function potential. 

We can calculate the leading-order influence of a finite- 
range potential by including the energy- dependent zero- 
range pseudopotential of Eq. For the TV-boson 



o 
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_ (a) N=3 


T V 
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- (b)N=4 
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0.01 
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FIG. 5: (Color online) Analysis of the three-boson [panel 
(a)] and four-boson [panel (b)] energies including scattering 
length and effective range effects. All energies are scaled 
by [a((jj)/a t (0)] 2 to emphasize the corrections due to finite 
range effects. The numerically determined finite-range ener- 
gies Efk are calculated for the Gaussian potential with spa- 
tial width ro = 0.005ct(cj) (red circles) and ro = 0.010cr(u;) 
(blue squares) , respectively. The dashed line shows the scaled 
perturbation theory prediction E^ + E^ for a zero-range, 
delta- function potential. The divergence at a t (0)/a(uj) = is 
due to the divergence of the effective range at zero scattering 
length. The unsealed energy shift vanishes when at(0) = 0. 
The solid lines show the scaled energies E {2) + E (3) + £ (1 ' 2) , 
which include the perturbatively calculated finite-range cor- 
rection £ (1 ' 2) . 



ground state, the pseudopotentialgives to first order in 
g' 2 an energy shift E^ 1 ^ [see Eq. Q]. At this order, the 
addition of V 2 is equivalent to replacing a t (0) by af(£ , re i), 
with Eq. ( [18] ) evaluated at the relative zero-point energy 
E re \ = 3/2jm units of Huj) of two non- interacting bosons 
in the trap. 

The solid lines in Fig. [5] show (E^ + £ (3) + 
E^)[a(u)/a t (0)} 2 as a function of a t (0)/a(u) for TV = 
3 and N = 4 trapped bosons, respectively. Combin- 
ing the perturbative predictions for zero-range contri- 
butions E^ + E^ and the effective-range correction 
gives excellent agreement with the nonperturba- 
tive finite-range energies. The comparison validates the 
perturbation theory and predictions derived in this pa- 
per for effective interactions including finite-range correc- 
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tions, through third order in a t (0). It also shows that the 
divergences in Fig. [5] at a t (0) = are due to the diver- 



gence of the effective range shown in the inset of Fig. Ill 



Finally, we note that the energy shift is proportional to 
the volume r e ff[a t (0)] 2 and goes to zero at a t (0) = 0, as 
expected. 



IV. FIRST- AND SECOND-ORDER EFFECTIVE 
INTERACTIONS 

A. Hamiltonian and renormalization condition 



The numerical results in Sec. |III| show that finite-range 
effects are important at third order in perturbation the- 
ory for realistic bosons. We incorporate these corrections 
by modeling the pairwise collisions of ultracold bosons by 
combining the zero-range pseudopotential 

V 2 (ri - r 2 ) = 4tt— a bare 5 (3) (ri - r 2 ), (20) 

where abare is now identified as the bare scattering length, 
and the effective-range potential 



V2O1 -r 2 ) = -^ 2 ,bare 

x fe (3) (ri - r 2 ) + *( 3 >(n - r 2 )^? 2 ], (21) 
which has the bare coupling constant 



#2,ba 



(22) 



The interactions of N ultracold neutral bosons can be 
described in quantum field theory with the Hamiltonian 
T~L = where Ho is the single-particle Hamiltonian 

and 



Hi 



l - J ^( ri )^t( r2 )[v 2 ( ri -r 2 ) 

+ V£(ri - r 2 )]^(ri)^(r 2 )dndr 2 . (23) 



The field operators ^(r) and ^(r) respectively annihilate 
and create a boson at position r. We assume the absence 
of intrinsic three- or higher-body interactions. 

The bosonic field is expanded over isotropic harmonic 
oscillator states with frequency uo as 

^(r) = ^ <l>nlm(r)a>nlm = ^ ^W^' ( 24 ) 
nlm i 

with hi annihilating a boson in orbital 0i(r). In the 
following we use the shorthand notation i = {nlm}, 
denoting the (dimensionless) single-particle energies as 
Si = £nim = (2n + I + 3/2), where n,l = 0,1,2,..., 
and i — is the {nlm} = {000} single-particle vibra- 
tional ground state. Substituting Eq. (24) into H and 



dividing by Huj, we define the dimensionless Hamiltonian 
H = H + Hi + H{, where H = £\ Sia\a u 

H\=\ \^zr^) ^K ij]k ia\a]a k ai, (25) 



2 \a(uj) 



ijkl 



and 



H i=2 (2 a(cu)3 ) Z^ K ^M a Wj a ^i' (26) 



The matrix elements 



K ij;kl = 47t[<jH] 3 J 4>*(r)4>*(r)Mr)Mr)dv 



(27) 



and 



^[a{uj)f [ [^(r)^(r)]V^fe(r)Mr)]rfr 



(28) 

are normalized such that i^oo ; oo = VV 71 " anc ^ ^oo-oo = 
(3/4) y/2/ir, with the semi-colon separating initial and fi- 
nal states and V? = (^? + ^?)/2. The factors of [a(cj)] 3 
and [cr(c<j)] 5 make the matrix elements dimensionless and 
lj- independent. As explained in Sec. [Hj we assume a 
regime where H[ can be treated as third order in per- 
turbation theory. 

The noninteracting ground state containing N bosons 
in the i = (i.e., nlm = 000) vibrational ground state 



is 



I TV) 



\0)/Vn\, with energy = Ne and 



So = 3/2. First-order perturbation theory in Hi gives 
eW(lj) = (1/2)N(N - 1)/7 2 (1) with 



(i) 



a. 



(l) / «b, 



(29) 



using (N\ a\a\aoao \N) = N(N — 1) and recalling that 
\N) denotes N bosons in the non-interacting vibrational 
state nlm = 000. The two-body, first-order coefficient is 



At higher orders in Hi, there are divergences due to 
the S- function potential (see e.g. [671 EH] ) - We regulate 
these by either truncating sums over intermediate states 
at a high-energy cutoff hu c , or by using an exponential 
regulator function. The former is more convenient for 
numerical approximations, while the latter is more con- 
venient for analytic results. In either case, we find at 
(2) 

second order that Z7 2 diverges as v /cJ^ and renormal- 
ization is required. Although this can be done using 
bare perturbation theory, in which infinities are absorbed 
by appropriately redefining bare parameters, we use the 
method of renormalized perturbation theory which pro- 
vides a systematic and self-consistent approach for calcu- 
lations beyond second order involving multiple divergent 
terms. 

Renormalized perturbation theory (e.g., see [33]) re- 
expresses the bare scattering length as 



«ba 



a t (uo) + a ct (u) ). 



(30) 
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Relative length scales 

a(oo ) = (/z/m A co ) 12 h.o. length at cd 
a(co) = (fi/rn A (ti) m h.o. length at co 

' a t (od )= (/z/m A co s ) 1/2 trap scattering length at co 

a t (0) = a^O) free-space scatt. length 

effective range 
width of Gaussian potential 



'eff 



\ 



{ l c = (fi/m A (D c ) 




1/2 



cutoff length scale 



FIG. 6: The length scales of interacting, harmonically 
trapped, ultracold bosons. We assume a separation of length 
scales l c <C at(cJo) <C cr (cjo), or equivalently a separation of 
energy scales uj c ^> cj s ^> cjo, where uj s — h/mAat(wo) 2 . The 
harmonic oscillator lengths ct(cj) and a (cjo) are of the same or- 
der, although a(uJo) is not necessarily larger than <j(uj). Sim- 
ilarly, we assume that a t (c<;o), a t (0), r e ff, and the Gaussian 
width ro are of the same order. The order of length scales 
within a group is arbitrary. 



A renormalization condition defines a t (uo) as the physi- 
cal scattering length for two bosons in a trap at frequency 
ujo. The cutoff dependent remainder a ct (co>o) is called a 
counterterm. For brevity, this notation suppresses the 
dependence of a ct (uJo) on uj c . In the following, we call 
a t (uJo) the "trap scattering length" at frequency co>o, to 
distinguish it from the energy-dependent free-space scat- 
tering length af(E re \) defined in Eq. (17). In the limits 



of zero relative collision energy and ujo = 0, the trap and 
free-space scattering lengths are equal, i.e., a t (0) = af(0). 
With the combination of V2 and V 2l the trap scatter- 
ing length a t (o;o) includes both the effects of the ujq- 
dependent dressing by quantum fluctuations to higher 
orbit als and finite-range effects. Note that the trap scat- 
tering length a t (co>o) does not, in general, equal the free- 
space scattering length a^(E re \) defined in Eq. (18) be- 
cause the latter does not correctly capture the influence 
of the harmonic confinement on the quantum fluctuations 
to higher orbit als. 

Together with the renormalization condition, the other 
key ingredient in renormalized perturbation theory is 
that the leading-order scattering length counterterm 
a ct (uJo) is proportional to [a t (uJo)] 2 ; in other words, it 
is a second- and higher-order contribution. This, plus 
the renormalization condition, systematically reorganizes 
the perturbation theory, order- by-order, so that it is an 
expansion in the physical value a t (uoo) instead of abare- 
Figure [6] summarizes the relationship between the char- 
acteristic length and energy scales for our model system 
of trapped ultracold bosons. 



Substituting Eq. ([30]) into Eqs. ([25]) and ([26]) gives 
Hi(uj;uj ) = V(uj]ujq) + V'(uj]ujq) + V c t(uj;uj ), (31) 
where the zero-range and counterterm operators are 

V(uj]uj ) = ^ \ ^r\- ) ^Kij. k ia\a\a k a h (32) 
Vctfauo) = \ ( J ^Kij. k ia\a]a k ai, (33) 

1 v j ijki 



and the effective-range operator is 
1/1 r e ff[a t (a;o)] 2 



2 \2 [(t(uj)} 3 



a k a t 



ijkl 



o 



r eS [a t (0)f 



(34) 



The renormalized perturbation theory is then organized 
based on the observation that V(uj]ujq) is proportional 
to a t (cJo), V c t(oj]UJo) is proportional to [a t (co>o)] 2 , and 
V'(uj\ ujq) is (for the regime considered here) proportional 
to [a t (uJo)] 3 . The single counterterm operator V ct (uj;uJo) 
cancels all divergences from the operator V(uj]ujo), at all 
orders in perturbation theory. In contrast, the effective- 
range operator V'(lj]LJq) leads to a nonrenormalizable 
field theory with the consequence that new counterterm 
operators are required at every order in perturbation the- 
ory beyond first order in g' 2 \ because we are only working 
to first order in g 2 in this paper, no additional countert- 
erms are needed. 

Note that the frequency ujq at which a t (co>o) is defined 
and the trap frequency uj for which we want to compute 
energies are independent. In the overview, we summa- 
rized our results for the special case where ujq = 0. The 
general case of arbitrary ujo facilitates renormalization of 
the perturbation theory. More importantly, the renor- 
malized perturbation theory is "calibrated" to a mea- 
sured value of at at a desired trap frequency co>o, and 
is then used to predict energies for trap frequencies uj not 
generally equal to ujq. 

We can now compute the ground-state energy 

E(uj; uj ) = e N + ^U 2 (uj; uj )N(N - 1) (35) 
+ ^U 3 (cd-cdo)N(N-l)(N-2) 
+ ^U 4 (uj; uj )N(N - 1)(N - 2)(N - 3) + .... 



We have used the semi-colon notation in Eqs. (31), (32), 
(33), (34), and (35) to distinguish between the roles of 



the frequencies uj and ujq. Before renormalization, the 
interaction energies U m (uj] ujo), found from perturbation 
theory in Hi(uj;ujq) : are functions of a t (uoo) and a ct (uJo). 
The renormalization condition can be expressed as 



U 2 (uj = ujo;ujq) 



flt(^o) 
a(uj ujq) 



(36) 
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which, in practice, is solved for a ct (cJo) to the desired 
order in perturbation theory. Another way of describ- 
ing the renormalization condition is that a ct (co>o) is tuned 
such that the first-order result is exact and the second- 
and higher-order corrections to the two-body energy van- 
ish when evaluated for two bosons in a trap with uj = uo . 
After renormalization, the interaction energies U m (uj; ujo) 
depend only on a t (uJo) and, moreover, the ^-dependence 
of the ground-state energy satisfies E(uo]uoq) = E(uj;u f ), 
for any pair of frequencies ujq and uj f . 



B. Energy at first-order in scattering length 

We use renormalized Rayleigh-Schrodinger (RS) per- 
turbation theory to compute the iV-boson ground-state 
energy E = X^ n =o E^ n \ where E^ is proportional 
to [a t (cJo)] n - We separate the contributions at each 
order into m-body energies, such that U m (uj;uJo) = 
^2 n Um\uj; ljq). The zeroth-order term is E^°\u) = 
eoN. The first-order energy shift is 

EW(w;w ) = (JV|Hi(w;wo)|JV) = (N\V(w;u> )\N) 

- 1 [E ( a t( w o)' 

~ 2 V 7T 



N(N-l), 



(37) 



using the fact that V, F ct , and V are 0(a\ i (uj^) j g(uj)\ 
O([a t (ojo)/cr(uj)} 2 ), and 0([a t (oJo)/cr(uj)} 3 ), respectively, 
and (N\ alala a \N) = N(N - 1). 

Comparing to Eq. (35), we see that the two-body en- 
ergy to first-order for any uj and ujq is 



with 



U^\uJ]UJq) 



,(1) 



S 1 ) 



(38) 



(39) 



For a trap with uj = cjo, the renormalization condition 
says that U^i^o] ojq) = y/2pn [a t (uj 'o) / & (ujq)] is the ex- 
act two-body energy. For uj ^ co>o, 11^(00 ^ ujq^ujq) is 
the leading order contribution to the full two-body en- 
ergy Z72(cj;cjo), but, as shown in the following sections, 
there are higher-order corrections that become increas- 
ingly important the more uj differs from ujq. 



ground state, Asij = S{ +ej — 2e , Zij is a normalization 
factor, and ij ^ 00 denotes summing over all z, j except 
z = j = 0. Equation (40) is modified from the usual 



RS perturbation theory because of the presence of the 
O([a t (ujo)] 2 ) interaction term V ct: which generates the 
counterterm contribution. 

The sums over intermediate states \ij) exclude the 
ground state z = j = 0, and are regularized using ei- 
ther a hard cutoff Asij < uJc/uj, or an exponential reg- 



-13 

ulator Ae" 1 — » e - A ^^/^ c Ae"- 1 , where Huj c is a high- 
energy cutoff. In the limit uj c /uj 
are equivalent. 



oo, these regulators 



Using Eqs. (32) and (33), we have 



(41) 



\ v / / ij^oo,ki 13 

and the expectation value is with respect to the non- 
interacting ground state \N — 2) oc aodo\N). The nota- 
tion ij 00, kl indicates that the sum is over all z, j, /c, I 
except z = j = 0. Wick's theorem gives 



(aid, Jala]) = 4(:a i a j a^aJ:) + 2(:a i a j a\a\:) 

= 45 ik (N-2) + 25 ik 5 jh (42) 

where :: denotes normal ordering, uncontracted indices 



are set to zero. 



t 

and contractions = 5^. Also, we 



have used (la^a+aj:) = 0, (:a\ M :) = (N - 2)(N - 
3)...(N — M + 1), and combined equivalent terms. Be- 
cause of the factor N — 2, the first term of Eq. (42) can 



be understood as leading to an effective three-body in- 
teraction, whereas the second term is a correction to the 
two-body interaction. 

(2) 

The second-order interaction energies (uj;ujq) and 



(2) 

C/g '{lo;lo ) can be extracted by evaluating Eq. (41) and 



comparing with Eq. (35). This gives 



X <> ( 43 ) 



C. Energy at second-order in scattering length 



The second-order energy shift is given by 

£ (2) = y ct;00 , 00 - E y °°f^ 00 , (40) 



zj^OO 



where V^/ = (zj|V|/c/) and Kt ; ij,/d = (zj|Kt|^)- Tne 
notation \ij) = Z^-aJaJaoao 1^) denotes the state with ei- 
ther one or two particles excited from the non-interacting 



and 



U^\uj;uj ) 



( 2) / otw y 
3 UhJ 



(44) 



The expressions for /3^\uj) and are defined in Ta- 
ble |TI| which also shows the explicit values calculated 
in Appendices [A] and [B] for an isotropic harmonic trap. 
We use the notation that and are associ- 

ated with n th -order, m-body processes. The sum that 
gives f3^ (uj) diverges with cutoff as ^uj c /uj, where uj c /uj 
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is approximately the number of harmonic oscillator lev- 
els included in the sum as a function of cj, for a fixed 

(2) 

cutoff uj c . The coefficient a 3 is convergent and in the 
limit uj c /uj — » oo is independent of uj. In the following, 
we only indicate the explicit uj dependence for coefficients 
that remain sensitive to uj in the limit uj c /uj — » oo, e.g., 

(2) (2) 

we write (uj) but a 3 . We use a hard cutoff to nu- 
merically evaluate the coefficients , af\ and af\ (see 
Sec. [V] for the definitions of the third-order coefficients.) 
For the coefficients a^\oL±\, and ol§ \ we find analytic 
results in the limit uj c /uj —> oo. Finally, using the ex- 
ponential regulator, we obtain analytic results for the 
coefficients fffi (uj) , f3^\uj) , and /3^ 3 \uj) for any uj c /uj. 
Equations (|43|) and (|44|) have also been represented 



diagrammatically, with factors of Kij. : ki[at(uJo) /a(uj)] as- 

signed vertices *X*? and contractions a^a^ representing 
excited particles assigned dashed lines Uncontracted 
operators a$ (or a\) are assigned incoming — • (or out- 
going • — ) lines. The count erterm is represented as 
Kij;ki[act(uo) / cr(uj)] = Intermediate states have one 
or more excited particles and contribute an energy de- 
nominator 1/Asij. For example, the diagram ^^<; is a 

graphical representation for the term [a t (uJo)/a(uj)} 2 , 



and U^ 2 \uj;ujq) 



6^"^<;. We obtain combinatorial 



(2) 

prefactors [e.g. —6 for ; (uj : ujq)] from Wick's theorem 
by counting the number of equivalent contractions, di- 
viding by 2 for every factor of a t (uJo) or a ct (uJo), and 
multiplying by ml for an m-body term. 

The renormalization condition through second order is 
U 2 (uj;ujo) = U^(uj -uj ) + u¥\uj ;ujo) + O([a t (uj )} 3 ) + 

O(r e fi[a t (uj )} 2 ) = /7 2 (1) (cj ;^o), and hence U^\uj ; ujo) = 
0. Diagrammatically 



XI 



xx 



Solving for the counterterm gives 

a c t(wo) = — I — — - I o 



(45) 



(too). (46) 



Substituting into Eq. (43) gives 



U^\uj;uj ) = cf\uj,uJo) 



Qt(^o) 
(j(uj) 



where the function 

4 2) (cj,cj ) = y/uo/uffiiuo) - ^ 2) (uj) 



(47) 



(48) 



can be used for any uj. (We have used a(uJo)/a(uj) = 
^uj/ujo above to simplify the expressions.) 

The form of the expression for the coefficient 

(2) 

c 2 (uj^ujq) ensures that the divergent terms cancel. For 
an isotropic harmonic oscillator, we show in App. |B 6[ 
using an exponential regulator, that 

$\w) = (2/tt) [ V W2^-(l-log2)]4-0(l/w c 1 / 2 ), (49) 



and thus 



,(2) 



(w,w ) = (2/tt) (1 - log 2) [l - v4^7 



(50) 



The renormalization condition is automatically satisfied 
since ^(ujq^ujq) = 0. For the special case when ujq = 0, 
we find 



(uj,0) 



,(2) 



(2/tt) (1 - log 2) = 0.19535.... (51) 

For brevity, we define c$ without arguments as the co- 
efficients Cm\u,0) for the special case when ujq = 0. In 
this limit, the coefficients are independent of uj. 

Combining the first- and second-order contributions 
for the two-body interaction energy gives 



U 2 (uj;ujq) 



(1) ( a t (uJo)\ ( 2 ) 



O 



a(uj) 
at(^o)] 3 
[a(uj)}* 



} (uj,uj ) 



a(u>) 
r eff [a t (wo)] 2 



O 



r{oj)f 



(52) 



(2) 

The coefficient a\ in Eq. (|44|) is finite and does not re- 



quire a regulator. For the three-body interaction energy 
we obtain 



U 3 (uj]ujq) = c ; 



Qt(^o) V 



O 



(2) 

r e fi[a t (uj )} 2 



O 



(53) 



where c 



(2) 



-6a 



(2) 



[a(uj)]t 

-0.85576... This value was previ- 



ously obtained in [18 , and is also calculated in App. B 1 



V. EFFECTIVE INTERACTIONS THROUGH 
THIRD ORDER 

We now extend our analysis to third order in the scat- 
tering length a t (uJo)- This is necessary to obtain the 
leading-order effective four-body interaction. Including 
the counterterm and effective-range interaction, the for- 
mula for the third-order energy shift is 



c/u 



Voo,ijVij,kiVki,< 



£( 3 )( w;wo )= rw ' l:,y,: ' 



00 



AsijAski 



(54) 



Vr 



00,00 



w c /w 



^oo,zj^i,oo 



U c /U} 



Ae 2 - 



ct;00,r/'^j,00 T/ / 

r ^00,00 • 



ASij 



The first term on the right-hand-side of Eq. ( 54 ) gives 
<t(w) ) 



£ 

ij^00,klqr,st^00 



KoO;ijKkl;qrKst;00 t 

Ae^Ae st 



(55) 
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FIG. 7: Sequence of boson-boson interaction induced transitions to higher orbitals. This example generates corrections to the 
ground state energy that can be viewed as an effective three-body interaction. The process, which involves three interaction 
vertices, arises at third order in perturbation theory and gives the energy shift [a t (coo) / cr(uo)] 3 derived in the text. Links 
labelled Vij represent intrinsic 2-body interactions between particles i and j. Black arrows represent virtual transitions to 
and from excited orbitals. Solid and dashed lines represent atoms in ground and excited vibrational states, respectively. The 
diagram on the far right shows the perturbation theory diagram for this process. 



where the expectation value is with respect to the non- 
interacting ground state with (N — 2) bosons. Applying 
Wick's theorem, Eq. (55) expands as 



ot(^o) 



1 

3! 
1 

4! 



Q/3?\u)N(N-l) 



(56) 



[12c4 3) + 12^ 3 \lj)}N(N - 1)(N - 2) 



[48a ( *l + 48a^ + 6a { 4 %]N(N - 1)(N - 2)(N - 3) 



,(3) 



.(3) 



+ ^60a {3) N(N - 1)(N - 2)(N - 3)(N - 4)) . 

We find effective two-, three-, and four-body interac- 
tions from the terms with four, three, and two con- 
tractions, respectively. The zero-contraction term van- 
ishes since (:aiaja\a\ai'aj'a\^\,:) = 0. Comparing to 
Eq. ( |35| ), we see that there is a two-body contribution 

(uj) [a t (ljq ) / a (uj)] 3 = VjtY, there are two three- 
body contributions 12c4 3 ^ [a t (uJo) /a(uj)} 3 = 12 *v** and 
12(3^ 3 \uj)[a t (uJo)/a(uj)} 3 = 12> :: *x, and so on. The def- 
initions for the coefficients /3r^\ /3^\ etc., are given 
in Table |TTJ along with the associated diagrams, asymp- 
totic behavior, and explicit forms for an isotropic har- 
monic oscillator potential (calculated in Appendices [A] 
and [B]) . Figure [7] illustrates one of the sequences of vir- 
tual transitions giving rise to . 

We next use Wick's theorem to evaluate the second 



term on the right-hand-side of Eq. (54), finding 



1 J3) AT 2 



-a 



4,3 



n 2 (n- ly 



-a^N\N- 



1) 2 (N -2) 



Qt(^o) V 
v(u>) J 



(57) 



where af\ and ot^ already appear in Eq. (|56|). Equa- 



4,3 



tion ( |57| ) can be separated into m-body contributions 
by expansion into terms proportional to N (N — 1), 
N(N — l)(iV — 2), etc. 



It 

,(3) 



is surprising, at first sight, that af\ an d 
several effective multi-body ener- 



a§' contribute to 



gies. From Table |II| a^[a t (cj )/cr(^)] 3 = and 

Qg 3 ^ [a t (uJo)/a(uj)} 3 = look like processes requiring 

four and five distinct particles, respectively. They also 
appear to be composed of "disconnected" sub-diagrams. 
In RS perturbation theory, however, the second term in 
Eq. (54) can be reinterpreted in terms of particles going 
"backward" in time (right to left), or alternatively an in- 
terpretation can be given in terms of holes. For example, 
the term —a ( fl[a t (uJo)/a(uj)} 3 gives a two-body contribu- 
tion if we view the two particles first going forward in 
time (left to right), colliding to an excited intermediate 
state, colliding back to the ground state, and finally going 
backward in time and colliding a third time. Diagram- 
matically, this can be represented by the connected dia- 
gram ^^S>J. Similarly, if only one particle goes back in 
time, it can collide with a third particle, giving the three- 



-6 



body contribution — 6a^l[a t (uJo) /a(uj)} 3 = 
which is also connected. In this paper, these related two-, 
three-, and four-body diagrams have the same numerical 
value: <> = <: > = . 



The third term on the right-hand-side of Eq. ( 54 ) gives 



the two- and three-body counterterm contributions 



~$\w)N(N - 1) + ±2a™N(N - 1)(N - 2) 



or 



and 



Oct(wo)\ / «t(^o) \ 
a(co) J \ a(co) J 



/3?\co)[a t (uj )a ct (uj )/a(uj) 2 ] = »: 



a 3 2) [«t(w )a ct (w )/o'(w) 2 ] = ^>c 



(58) 



(59) 



(60) 



These counterterm contributions, shown in Table [ill can ~ 
eel the divergences from yCX^, ^OOC and x- The 



disconnected 5-body contribution from Eq. (57), gener 



ated by the term with a single contraction, cancels with 
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Energies (Diagrams) 


Coefficients 


Asymp. 


Isotropic H.O. coefficients (ou c — > oo) 


l st -order in £t = at(wo)/cr(uj) 


4% =X 


ol^ = Koooo 


N.A. 


J% = +0.797885... 


2 nd -order in f t 


a (2)e2 _ . 
"3 St — / 


„,( 2 ) _ \^ ^OOOi^iOOO 

3 - 2^ A£i0 


a + e Wc / W 


(|) + log(8 - 4v^3) - 1] = +0.142626... 


&\u)e t = >::< 


p {2) ^OOij^ijOO 
02-2. A £ij 


/ 

V w 


(Di^-a- 1 ^ 2 )-!^] 


3 rd -order in £t 


^\^e t = >;:*X 


„(3) ^ KQQijKijklKj^lQQ 


it) 




^V)*? = <:\ 


/q(3) K 00ij K ij0k K k000 
P 3 _ 2. A £ij A £fc0 


\j U> 




4 3) £? = >r 


(3) ^OOij^jOOfc^ifcOO 
«3 ~2^ Ae, 3 Ae tk 




+0.56494 ± 0.00001 (estimate) 


a (3).3 _ 0~C 


^(3) Kooij^jOOO^iOOO 


a + e -"c/" 


+0.077465... (numerical) 


a (3) £ 3 - 


(3) KoooiKiooj Kjooo 
4,2 - 2^ A £i0 A£ i0 


a + e" Wc / w 


+0.051099... (numerical) 




;~( 3 ) ^ooii^oooo^ijoo 
^4,3-2- A £ 2 . 




(l) 3/2 III + lo § 2 - 1 ( lQ g 2 ) 2 ] = +0-438946... 


4 3) «? = 


~,( 3 ) _ ^OOOi^OOOO^iOOO 


a + e -u c /uj 


4(27r 3 3/2 4^3 (1, 1, 1, 5/2; 2, 2, 2; 1/4) = +0.051916... 


Counterterms through third order 


- , 4 i} xct - >:;, /3< 2) (cu )X c t s t - >::<, a ^ Xctit - x< 


Leading-order effective range terms 


4 1>2) (3Sj) = X- 4 1>2) = ^ooo = ! ) V2 = +0.598413... 


Other relations: 


>c*C XX XX? 

# (four-body) = ^^^^ (three-body) = 


(two-body) = «i 3 ^ t 3 


>•<*< ^^S^c: > *><i 

^^♦^^ (five-body) = ^^.^^ (four-body) = 


(three-body) = <4 3) £j? 



TABLE II: The coefficients for all interaction processes contributing to the two-, three-, and four-body interaction energies 
through third-order in perturbation theory in £t = at(uJo)/cr(uj). The first column shows the diagrams from which the ra-body, 
n th -order coefficients affl and (3m\w) can be reconstructed. The coefficients as multidimensional sums are given in the second 
column. Sums are over all indices i, j, k, ... except combinations that give a zero energy term in the denominator. The third 
column gives the asymptotic behavior of the coefficients in terms of the cutoff uj c and a constant a. The last column gives the 
explicit values for the coefficients for an isotropic harmonic oscillator potential. These values are obtained in the Appendices. 
The table also shows the counterterm processes, the leading-order effective-range contribution, and other relations needed for 
the renormalized perturbation theory. 



the disconnected five-body term in Eq. (56), and there is 



no effective five-body interaction at third order. Finally, 
the last term on the right-hand-side of Eq. ( 54 ) gives the 
effective-range contribution 



Qt (^o) 



£ -) ( ^ 0) =i»r^-i»(^)(^ 

(61) 

from which we extract the effective-range two-body in- 
teraction energy 



y-y-(l,2) / \ 



(1,2) ( r e ff[a t (^o)]^ 



X ( 62 ) 



The coefficient a^ 1 ' 2 ^ is given in Table [III The special case 



gives Eq. (JT)) and Eq. ([8). 



A. Two-body interaction energy 

Adding all two-body contributions through third or- 
der, we obtain 

U 2 (uj; ujo) = X - X< + X - + >X< 

-2 O +X + ^(«t) (63) 

(i) /«ctM\ ( 3 ) ( at(u} ) \ 3 
, a,2) / reff \ / a t (oJo)\ 



(64) 
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Note that all diagrams in Eq. (64) are connected, when 



interpreted in terms of both forward- and backward prop- 
agating particles. All coefficients are given in Table [TT| 
For brevity, in Eq. (64) and the following we adopt 



the convention that 0[a t 4 ) means O([a t (ujo)/a(uj)] 4: ) 

O(r eff [a t ( W0 )] 3 /^(w) 4 )- 

The counterterm, found in the previous section to 
second order, must now be recalculated using the 
renormalization condition through third order. This 
adds a third-order term which cancels the divergence 
from f32 3 \(ju)[a t (uJo)/(j(uj)] 3 = YxV, as well as the 

(1 2) 

effective range contribution [/^ ' (uj;ujq). Solving the 
renormalization condition L^f^^o) = U^\ujq; ujq) + 



U^\ujq] c^o), and hence U^^uq] ujq) + [/^(^oi^o) 

(1 2) "~ 

U 2 ' j (uq;ujo) = 0, we find a ct (uo) from 



(2)/ 



zero-range interactions found in [34 . This agreement 
between the quantum mechanical and quantum field the- 
ory solutions is a nice illustration of how the renormal- 
ized effective field theory captures the correct low-energy 
physics. Interestingly, if r e ^ ^ 0, our result for [^(^O) 
still agrees with the solution in [34], if that solution is 
Taylor expanded in a t (0) and r e ff[a t (0)] 2 after making 
the substitution Of(0) — )> af(0) + (l/2)r e ff[af(0)] 2 /c 2 el , pro- 
viding further evidence of the universality of the higher- 
order perturbative results derived here. 

Another important special case is a; = uo. Since 

4 2) (cj ,^o) = 4 3) (^o,^o) = d£' 2 \(jJo,LJo) = 0, the pre- 
dicted two-body energy is 



U 2 (u> = <*>;<*>) = <g ) f^l) +0(4), (68) 



4 1} (^W>m(^Y 

V <K^o) / V ^(^o) / 



4 3) (^o)-2[4 2) (c )] 2 /4 



* (3) 

*4,3 



Qt(^o) V 
a(uj ) J 



7 (1>2) 



S)(^t) +0 «'- 



(65) 



Diagrammatically, this can be expressed as 



x - 2>X = >x + 



<:*:> X ( 66 ) 



with all diagrams evaluated at uj = cuo. By including 

(1 2) 

U 2 (uj;uo) in the counterterm equation, the renor- 
malization condition for a t (oJo) includes both zero-range 

and effective-range contributions. If we do not in- 

(1 2) 

elude U 2 ' (a;; cjo) in the renormalization condition, then 
a t (oJo) is the trap scattering length for zero-range poten- 
tials. Substituting the counterterm from Eq. (65) into 



Eq. M) and using ^ (u) = [^ 2) (cj)] 2 , which is 



proven in Appendix |B 7\ we find after some algebra that 



Qt(^o) 



TT ( \ (!) f a t(^o)\ . (2)/ \ 

+ ^' 2> <— >(^))(^r) 2+oM) - (67) 



reproducing the renormalization condition that a t (uJo) is 
the physical trap scattering length for two bosons at fre- 
quency UJQ. 



Frequency-Dependent Effective Interaction Coefficients 



: (2/7T) 



1/2 



Two-body 

(1)/ \ (1) 

4 2) (u, ljo) = (2/tt) (1 - log 2) [l - vWT^ 
4 3 V,a;o) = (2/tt) 3 / 2 (1 - log2) 2 [l - yWT^] 2 

- (2/tt) 3 / 2 (tt 2 /24 + log 2 - § log 2 2) [1 - wo/w] 
4 1,2) (w,wo) = 4 1,2) [1 - ^oH = (3/4) (2/tt) 1 / 2 [1 -cjoH 
Three-body 



4 \uj,ujq) = -6a 



(2) 



= -12a< 2 >c< 2 > (^ol/aW 
-[12c4 3) - 6a^ - 18a^ 3) ] 



Four-body 



(3) 



-72a| 



(3) 



TABLE III: The functions (a;, cjo), which determine the 
n th -order contributions to the ra-body effective interaction 
energies, and d^' 2 ^ (cj, cjo), which determines the leading-order 
effective-range correction, for neutral bosons in a harmonic 
potential of frequency cj, in terms of the scattering length 
at(cjo) defined at trap frequency ujq. The special case ujq — 
reduces to the results given in Table 1. 



where c^\lj,ljq) and cS/'^fk;, ljq) are given in Table 

Recall that in the formula U 2 {uj\ujq) the first argument 
u is the trap frequency for which we are interested in pre- 
dicting the two-body energy, and the second argument ujq 
is the trap frequency at which the two-body trap scatter- 
ing length a t (uJo) is defined or measured. The coefficients 
for uo = are given in Table [l| If r e ff = 0, these values 
reproduce through third order the exact solution for the 
ground state of two harmonically trapped bosons with 



7(1,2) 



III 



B. Three-body interaction energy 



In [18], we obtained the effective three-body interac- 
tion energy to second order. We now determine the 
next-order correction by combining all three-body con- 
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tributions through third order, giving 



U 3 (uj;oj ) = -6>C; 



12 



+ 12 



12 y<. 



«t(^o) 



18^^+0(4) 



_ 6a (2) y 12a ( 8) /".u-„p 



(69) 



+ i 2/ f > M f *£*T\ - 12af f f ^¥ 

(3) / a t (^o) V 1g (3) / «t(^o) \ 3 , 4n 



As anticipated, the two disconnected terms that depend 



(3) x:>c 

4,3 



(3) 

cancel and, at this order, c\ ' is inde- 



on a 

pendent of uj and ujo. The leading-order contribution to 
the four-body energy does not require renormalization, as 
is true for all leading-order m-body terms. Comparison 
of cf^ and c^p reveals, however, that they are of similar 
magnitude and therefore for a consistent and accurate 
treatment both corrections need to be included. Because 

(3) 

c 3 requires renormalization, we see why the systematic 
renormalization of divergences is needed even though the 
leading-order contribution to the four-body interaction 
energy could be obtained without these considerations. 



(3) 

Representing the three-body contributions from a\ 3 and 

a 5 using reversed (left to right) particle lines, as previ- 
ously described, we again find that only connected dia- 
grams contribute. 

For the three-body energy, it is sufficient to use the 



second-order counterterm in Eq. (46). In Appendix B8 

we show that @P (uj) = j^P {uj)ap /(dp. From these 
results it follows that the difference between the^individ- 
ually divergent contributions in Eq. (69), 12 ^*"*x anc ^ 



12 /^>c^, is finite. After some algebra, we find that 
tt ( \ ( 2 ) f a t(u )\ 2 ( 3 ) f a t (uj )Y 



+ OK 4 ), 

^p and cP ( 
If ujq equals zero, we find 



where c 3 and c 3 (uj,ujq) are given in Table 



III 



4 3) = -12(1- log 2)4 1) 4 2) 
= +2.7921 ±0.0001. 



[12a : 



(3) 
3 



6a 



(3) 
4,3 



(70) 



184 3) ] ©• In Sec 



(71) 



The error reflects a one standard deviation uncertainty 
due to the extrapolation of the n ume rical estimate for 

(3) 

a 3 to the limit uj ( 
case is uj = cjq, giving 



00 (see App. B 2). Another special 



12ai 3) - 6a 



(3) 
4,3 



18a 



(3) 



+3.2112 ±0.0001. 

(72) 



C. Four-body interaction energy 

Finally, we calculate the leading order contribution to 
the effective four-body interaction energy. We find 



-48^x < -72 ~r +0(at) 



with coefficient 



(73) 



cX> = 4Sa\J + 48a^ - 72a^ ; = +2.43317.... (74) 



VI. SUMMARY 

We have derived effective two-, three-, and four-body 
interaction energies for N bosons in an isotropic har- 
monic trap of frequency uj. These energies are functions 
of the trap scattering length a t (uJo) and harmonic oscil- 
lator length ct(uj)i and include both renormalization ef- 
fects due to quantum fluctuations to higher-orbitals and 
leading-order finite-range corrections. The frequency ujq 
at which the scattering length is defined plays a role 
closely analogous to the low-energy scale at which cou- 
pling constants are defined in high-energy effective field 
theories (e.g., see [33 j). The formulas for the interaction 
energies are given in Eqs. (67), (69), and ([73]), and are 



expressed in terms of the functions c^\lj](jJq) given in 
Table [Hi} In turn, these functions require the coefficients 



01$ and /3m\uj) given in TableJllJ The special case when 

UJ 



Jo = is summarized in Table [jand Eqs. (10), (11), and 



III we showed that these results give ex- 
cellent agreement to numerical simulations for ultracold 
bosons interacting through a Gaussian model potential. 

We find at third-order in a t (uJo) that the shifts to 
the effective three- and four-body interaction energies 
are comparable, showing that the renormalized three- 
body interaction needs to be taken into account when 
the leading-order four-body interactions are considered. 
In the future, we plan to use this formalism to deter- 
mine the effective multi-body interactions for other po- 
tentials, such as anisotropic traps or the anharmonic sites 
of an optical lattice. The cross-over from the perturbative 
small scattering length regime to the universal regime of 
Eflmov physics is also very interesting and diagrammatic 
resummation techniques can be used to study the on- 
set of nonperturbative behaviors. For example, collapse 
and revival experiments suggest that four- and higher- 
body interactions may be present in the data [17 , but 
our results also show that in these systems a t (0)/a(uJo) 
is large enough for significant nonperturbative effects to 
be important, and we would like to better understand 
this physics within our framework. A unified descrip- 
tion of elastic and inelastic interactions (e.g. three-body 
recombination physics [5j[57]) would also be useful. 
More immediately, the results in this paper can be ap- 
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plied to investigations of finite-range interactions, can be 
used for precision experiments probing for the possible 
existence of intrinsic three- and higher-body interactions, 
and can enable explorations of fundamental concepts 
in effective field theory including renormalization and 
energy-dependent (running) coupling constants. For ex- 
ample, the influence of intrinsic higher-body interactions 
would cause deviations from our predictions, which are 
based on only intrinsic two-body interactions. Moreover, 
we can engineer and exploit useful effective interactions 
using a combination of magnetic Feshbach resonances [2] 
and the ability to tune the few-body interactions by con- 
trolling the trap parameters and shape [53l [6Tt . One of 
our longer-term goals is to use this physics to develop 
nonlinear measurement techniques. For example, the 
nonlinear dynamics seen in collapse-and-revival experi- 
ments can lead to better than shot-noise measurements 
of the m-body interaction energies, or it may be possible 
to exploit strongly correlated non-equilibrium states in 
lattices for new types of sensing. In this way, the rich 
physics of renormalization and nonlinear quantum dy- 
namics could be used to create new types ultra-cold atom 
simulators, quantum information processors, or quantum 
sensors. 
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Appendix A: ^-function boson-boson interaction 
matrix elements for an isotropic harmonic trap 

This appendix derives interaction matrix elements for 
bosons in an isotropic harmonic oscillator trap with fre- 
quency uj and zero-range 5- function interactions. Alter- 
native methods for obtaining these matrix elements are 
given in [B21IZQ1. 



1. Isotropic harmonic oscillator wavefunctions 

The calculations are most conveniently performed in 
coordinates scaled by the harmonic oscillator length 
(j(uu). In spherical coordinates, the normalized, dimen- 
sionless isotropic harmonic oscillator states \nlm) have 
wavefunctions nZm (r) = (r\nlm) = Xni 0) Yj m (0 , 4>) , 



where Yj m (0, <fi) are spherical harmonics. The radial func- 
tions are 

Xnl (r)=N nl r l e- r2 / 2 Lil+^(r 2 ), (Al) 
where L„ (r) are associated Laguerre polynomials, 



2r(ra + l) 



T(n + 1 + 3/2) 
are normalization constants, and 

T{n + l + Z/2) 



4^/2) (o) 



r(n + l)r(Z + 3/2)' 



(A2) 



(A3) 



The single-particle ground state is </>ooo( r ) = 
7r -3/4 g -r /2 R eca n that we use the shorthand no- 
tation i = {nlm} for states with vibrational quantum 
number n, angular momentum Z, and angular momen- 
tum projection quantum number m. The single-particle 
energies are S{ = e n im = 2n + / + 3/2. A complete 
set of (un-symmetrized) two-particle wavefunctions is 
\ij) — \n\l\m\, n^ii'm'i). For convenience, we define the 
(dimensionless) two-particle energy differences 

A<S^j = A.£ ni i irni ,77,2^2^2 = ^nilimi H~ ^ri2^2^2 2^000 

(A4) 
(A5) 



2m + 2n 2 + /i +Z 2 - 



2. Matrix elements in the single-particle basis 



The matrix elements Kij-ki defined in Eq. (27) cor- 



respond to transitions \kl) —> \ij) with two-boson basis 
functions \ij) and \kl) from the |niZirai, 712/2^2) basis. 
In this subsection, we evaluate the subset i^oo of these 
matrix elements given by 

imi,n2/2^2;000,000 
= 47r J ^mhm, ( r ) <t>n 2 l 2 m 2 ( r ) ^000 0) 0000 0) dr 
= *Zi,i 2 <W,-m2^s.p.( n l> n 2jil)> ( A 6) 

where S a 5 is the Kronecker-delta and 



^s.p.( n l> n 2,0 



:N ni iN n2 i 



x j ^+1/2) ( r 2) L (/+l/2) ( r 2) e -2r\2l + 2 d ^ (A?) 

The subscript "s.p." means single-particle basis, and we 
have used the orthonormality of the spherical harmonics. 

We next use the complex contour integral representa- 
tion [71] 



L (J+l/2) (r 2) 



1 

2iH 



z/{l-z) 



{1-z) 



dz, (A8) 
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with a clockwise contour circling the pole at z — 0. Sub- 
stituting and then integrating over r gives 



K s . p .(nun 2 ,l) = —N nil N n2l T(l + 3/2)x 
V 71 " 

1 / dzi 1 / dz2 1 



*r 1 - ^ ^ (2 - ^ - z 2 ) 



Z+3/2 * 



(A9) 



Applying the Cauchy residue theorem twice, first inte- 
grating counter-clockwise around the pole at z 2 = 0, and 
then around z\ = 0, and substituting in the expressions 
for the normalization constants N n \ gives 



Appendix B: Perturbation theory coefficients 
through third order 

In this appendix, we compute for neutral bosons in an 
isotropic harmonic potential the m-body, n th -order co- 
efficients ol$ and needed for the perturbation 
theory through third order. We first evaluate the coeffi- 
cients otp , , , af\ •> a< t\ •> anc ^ •> which are finite 
and cj-independent in the limit that uj c /uj — >• oo. Then 

(2) (3) (3) 

we evaluate the coefficients (cj),/^ (cj), and (3% (cj), 
which diverge as u c /u —> oo. 



^s.p.( n l> n 2,0 = A/ - X 



(A10) 



+ n2 + i + 3/2) 



Vr(ni + i)r(n 2 + i)r(m + 1 + 3/2)r(n 2 + z + 3/2) ' 

This expression also gives the matrix element for the 
transition \0i) —> \0k). 



3. Matrix elements in relative and cent er-of- mass 
particle basis 

It is simpler to compute some matrix elements by 
switching to a basis of states = \nlm, NLM), 

with normalized relative and center-of-mass wavefunc- 
tions </>nZm(r) < £iVLM(R<) defined in terms of coordinates 
r = (r 1 — r 2 )/ a/2 and R = (i^ + r 2 )/ V% and (dimension- 
less) two-particle energy differences 



Ae nt 



Im.NLM 



2n + I + 2N + L. 



(All) 



Working in the \nlm, NLM) basis and using the 
fact that the interactions conserve the center-of- 
mass motion, the matrix elements for the transitions 

1^0 ~^ fij) are Kfj-ki = Knlm,n'l'm'',NLM,N'L'M' = 
Kreifan^SiflSrnflSi'flSrn'flS^N'SLtL'SMtM'i wnere 



K Te \(n,nf) 



2 ^oo(0)^oo(0) 
^ |0ooo(O)| 2 



(A12) 



only depends on the principle quantum numbers for the 
relative motion. Below we use the fact that K re \(n,n') 
factors as 



K ve \(n,ri) 



yJ^K Tel (n,0)K Tel (ri,0), (A13) 



and 



K rel (n,0) 



/r(n + 3/2) 



7r 3 / 4 y r(n + i) 

Also, K re \(0,0) equals ^2/^. 



(A14) 



(2) 



1. Three-body, second-order coefficient a 

In the single-particle basis \nilimi,n 2 l 2 m 2 ) ? the con- 
tribution has the coefficient 



.(2) - 



-KoO;Oi-?Qo;00 



(Bl) 



where the sum J^i/o * s over an a U° we cl single-particle 
states excluding the ground state. Due to angular mo- 
mentum conservation only i = {nlm} with I = m = 
contribute, and Aeio — 2n. Evaluating the sum gives the 
analytic result 



(2) _ ^ ifs.p.(^0,0) 2 

2n 




(B2) 

log(8 - Ay/S) - 1 I = 0.142626.... 



The terms in the summand become smaller exponen- 
tially with 77, and converges to the asymptotic form 
a + 0(e _u;c/u; ). This behavior will be true for all "tree 
diagrams" which, like ^^<;, have no closed loops. 



(2) 

2. Three-body, third-order coefficient a 3 

Continuing to work in the single-particle basis, the con- 
tribution has the coefficient 



a 



(3) 



k:00 



i?y00,i/c/00 



AsijAsik 



(B3) 



with i — {niZirai}, j = {n 2 l 2 m 2 }, and k = {773/37773}. 
Due to angular momentum conservation, we have l\ = 
h = h a nd mi = —7772 = —7773. Using Eqs. ( |A6[ ), ( A10| ), 



and ( A5 ) , the coefficient is given by 



4 3) = £ (2Z 1 + l)x 
K s . P Xn u n 2 , /i)i^ s .p.( n 2, 77 3 , h)K s . P X n i, n 3, h) 

(2m + 277 2 + 2/i)(277i + 277 3 + 2/i) 



(B4) 
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FIG. 8: (Color online.) Plot of numerical approximations 
to the sums which give the coefficients (red circles) and 
a< 4,i (blue circles), with a hard cutoff uj c in the energy of 
the intermediate states. The data is plotted versus {uj/uJc) 1 ^ 2 
on the bottom axis (the top axis shows the corresponding 
value of uj c /uj). The black lines are the least-square fits to 
the expected asymptotic behavior a + ^(cj/cJc) 1 / 2 + c(cj/cj c ). 
The values for and a 
to the y-intercept (u. 



(3) 

^4 3 are obtained by extrapolating 
oo). To estimate the one-standard 

J 3 ) 



deviation uncertainty in a 3 , for which we do not have an 
analytic value, we use the difference between the extrapolated 
and analytic values of ot^l- 



numerical determination demanding. We obtain an esti- 
mate by fitting numerical approximations versus uj/uj c to 
the asymptotic form a + b(u /ojc) 1 / 2 + c(uj/uj c ), dropping 
terms that are G[(uj/uj c ) 3 ^ 2 ]. The best-fit constants a, 6, 
and c give the curve a^\uo c ) shown in Fig. tal The best 
estimate for a K 3 , found by extrapolating uj c /uo — > oo, is 



a 



(3) 



0.56494 ±0.00001. 



(B5) 



To determine the one-standard deviation uncertainty in 

(3) 

«3 ; associated with our extrapolation method, we have 



(3) 

compared the analytic value of a\ 3 given in Eq. (Bll) 



(3) 

to the value for a\ 3 found by numerical extrapolation. 
The comparison is shown in Fig. [8] 



(3) (3) 

3. Four-body, third-order coefficients a\{ and a\ J 2 



The contributions and >$< give the coeffi- 

cients 



(3) _ K00;ijKj0;00Ki0;00 



and 



(B6) 



where the sum is over < 2ni + 2n2 + 2Zi < co> c /cj and 
< 2ni + 2ns + 2Zi < uj c /uo. The factor (2Zi + 1) arises 
due to the sum over the quantum number m\. 

f3) 

We have not fou nd a n analytic expression for a 3 , and 
the sums in Eq. (B4) converge slowly, making precise 



(3) _ ^ ^00;0i^i0;0j^j0;00 



M,2 



Asi As j0 



(B7) 



respectively. Due to angular momentum conserv ation , 
only z, j with / = m = contribute. Using Eqs. (A6), 
(A10), and (A5), we obtain the numerical results 



a W y g,.p.(n 1 ,0,0)^ p (n 2 ,0,0)^ p ,(n 1 ,n a ,0) = Q Q77465 (Bg) 

and 

a (3) = y g,.p.(ni,0,0)g..p.(ni,n a ,0)^ p .(n 2 ,0,0) =a051099 ^ 

These are tree-diagram processes, which, like a 3 , converge quickly, thereby making it is easy to obtain a precise 
numerical approximation from a small number of excited orbitals. 



4. Four-body, third-order coefficient a 



(3) 



The two-, three-, and four-body contributions < *; > . >i^> , and have the same coefficient, 



v (3) _ ^00;ij^00;00^i;00 _ ^00;fo^00;00^Qj;00 (Rift) 

^ As 2 . ~ ^ As 2 -. ' 1 ] 



where in the last expression, rather than evaluating ^000 ( r i ) </>ooo ( r 2 ) = </>ooo( r ) ^000 (R) and sum over relative 
the sums in the single-particle basis, we observe that 
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and center-of-mass bases states = \nlm, NLM), ex- 
cluding ij = 00. The interactions conserve the center-of- 
mass motion, implying j = {NLM} — {000}. Angular 
momentum conservation gives i = {n00}. Finally, using 

A£ n Q0,000 =_2n and i^0;00 = i^ n 00,000;000,000 = ^rel(^O) 



from Eq. (A14), we obtain the analytic result 



(3) 
*4,3 



2 [K re l(n,0)] 2 

n>0 



4n 2 



3/2 ^ 
^24 



log2 



^(log2) 2 ] 



(Bll) 



0.43894.... 



If we include the exponential regulator, we confirm that 



* (3) 



converges as (uj/uJc) 1 ^ 2 . Because the sums for a 



(3) 
4,3 



(3) 

and «3 have the sa me as ymptotic behaviors, we use the 



exact result in Eq. (Bll) to determine the accuracy of 
the extrapolation for shown in Fig. [s] 



(3) 



5. Five-body, third-order coefficient a 

The three-, four-, and five-body contributions 
and ^J^^ have the same coefficient 



(3) 

V v J 



E 



^00;Oi^Q)0;OO^Qo;< 



00 



(B12) 



Working in the single-particle basis |niZirai, ^2/2^2), we 
obtain the analytic result 



v (3) 



2 [K s . p .(n,0)] 2 



E 

n>0 



4n 2 



4(2tt) 3 / 2 
3/2 



4 F 3 (l,l,l,5/2;2,2,2;l/4) 



/2\ ' 1 

f-J [-Li 2 (l/2->/3/4)-log(l + >/3/2) 



i(log(l + v / 3/2)-log2) 2 + log2] 



0.051916..., 



(B13) 



where p F g is a generalized hypergeometric function, and 
hi2(z) is the polylogarithm function. Evaluation with a 
regulator function shows that this expression converges 
as (cj/c^c) 1 / 2 . 



(2) 



6. Two-body, second-order coefficient p£ 

The coefficients (uS) diverge when uj c /uj —> 00. The 
two-body contribution yC>C nas the coefficient 



4 2) M = E 

ijVOO 



w c /u; 

E 

2^0 



^00; £0^0; 00 

A^ 



^ ^o 

(B14) 

where we have switched to the relative and center-of- 
mass basis = \nlm, NLM) in the last expression. 
Using the fact that only I = m = and j = states 
contribute greatly simplifies the evaluation of P\ i^) by 
reducing the multidimensional sum to a single summa- 
tion. Using if n oo,ooo;000,ooo = K Te \(n,Q) from Eq. ( |A14[ ) 
and the exponential regulator l/Aem — > e~ 2n ^^ c) /2n, 
we obtain 



# 9) h = E 



n>0 



[i^rei(n,0)] 2 
2n 



2a; V 8 y 2 V 2w r 



(B15) 

+ o(iK). 



This coefficient diverges as ^uj c /uj, but as shown in the 
main body of this paper, the divergence cancels after 
renormalization, leaving a finite correction proportional 
to (2/tt)(1 -log 2). 



(3) 



7. Two-body, third-order coefficient /3. 

Next we consider the contribution yOO\, with coeffi- 
cient 



n(S) / \ _ KoO]ijKij]klKkl;00 ^Q0;?Q^0;fc0^fe0;00 (Rlfi") 

zj^00,/c^00 J ^o,/c^O KU 
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where we again switch to relative and center-of- mass basis states an d use the selection rules. Inserting exponential 
regulators for both energy denominators in Eq. ( B16| ) and using Eq. (A13), it follows that 



/o(3)/ \ K re i(n,0)K re i(n,n )K Te \(n , 0) _2(n / +») 

Po (uj) = > ; e 

2 v ^ 4nn' 



n>0,n'>0 




n>0 



[^rel(n 7 Q)] S 

2n 




n'>0 



In' 



(B17) 



This factorization result is important for the renormalizatio n of the two-body interaction at third- and higher-orders. 



8. Three-body, third-order coefficient fi^ (u) 



The contribution /*"*x gi yes the coefficient 



E 



KQ0;ij Kij;0kKk0;Q0 



E 



Kqo; m Km ; o/c i^fco ; op 
A^o^^/co 



(B18) 



In the last equality, we replaced the sum over single- 
particle intermediate states with a sum over rela- 
tive and center-of-mass states and then used the 
selection rule j = 0. The sum over k remains over the 
single-particle basis. We therefore require the "mixed- 
basis" matrix elements i^0;0fc- Using the selection rules 
I = m = for the relative motion, and l\ = m± = 
for the single-particle motion, we need only i^0;0fc = 
^mixed(^ ^i) with l = {n00} and k = {ni00}, where 



#mixed(n,rii) = ]J^(2tt) 3/2 J 4>* n00 (r) 05 oo (R)(j( 3 )(r)0 niO o(ri)0ooo(r2)drdR 

= y|(2^) 3 / 2 nOO (0) | 0* oo (R)0 niOO ( A)0 ooo (^) dR? (B19) 
and ri 5 2 = (R ± r)/y/2. Substituting in harmonic oscillator wavefunctions gives 

poo 

^mixed(n, m) = lGV^nOO (0) N ni0 / l£™ e - 2x \ 2 dx, (B20) 

./0 

where x = |R| />/2 and we h ave integrated over the angles. Noting that the remaining integral over x is proportional 
to if s . P .(ni,0,0) in Eq. ( |A7| ), we find that 

^mixed(n,ni) = y|K re i(n,0)K s . p .(ni,0,0). (B21) 

Inserting exponential regulators for each energy denominator in Eq. ( |B18 ), we obtain 

/o(3)/ \ ^rei(^,0)K m i xe d(n,ni)iC s .p.(ni,0,0) _2(n+n 1 )u; 

Po [UJ) = > e 

3 v ' ^ Ann x 

n>0,ni>0 

= V /|(E [W ^° ,0)ia ^ ) (E 1 ^*^-) =«rf. (B22) 

The factorization of /3g in the finite part /a^ and the divergent part fi\ (u) is important for the renormalization 
of the three-body interaction at third order. 
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